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Finite difference discretization schemes preserving a subgroup of the maximal Lie invariance 
group of the one-dimensional linear heat equation are determined. These invariant schemes 
are constructed using the invariantization procedure for non-invariant schemes of the heat 
equation in computational coordinates. We propose a new methodology for handling moving 
discretization grids which are generally indispensable for invariant numerical schemes. The 
O | idea is to use the invariant grid equation, which determines the locations of the grid point 

at the next time level only for a single integration step and then to project the obtained 
solution to the regular grid using invariant interpolation schemes. This guarantees that the 
scheme is invariant and allows one to work on the simpler stationary grids. The discretization 
errors of the invariant schemes are established and their convergence rates are estimated. 
Numerical tests are carried out to shed some light on the numerical properties of invariant 
discretization schemes using the proposed evolution-re-mapping strategy. 
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1 Introduction 

a 

Discretization schemes for differential equations that are not solely constructed for the sake of 
reducing the local discretization error as much as possible, but rather to preserve some of the 
intrinsic properties of these differential equations have become increasingly popular over the last 
decades. While preserving one of these properties, namely conservation laws, led to the design 
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£N) of conservative discretization schemes which are quite popular in the scientific community [4, 

18, 23] (and in particular in the geosciences, e.g. [17, 30]), there are other geometric features 
of differential equations that can be attempted to be preserved as well that have received less 
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attention from the side of numerical analysis so far. One of these features are symmetries of 
differential equations. While there have been theoretical advancements on the methodologies 
of finding numerical schemes that preserve the maximal Lie invariance groups of systems of 
J> differential equations over the past 20 years or so [12, 22, 29, 34], little is known about the 

numerical properties of these invariant schemes. A part of the problem is, while conservation 
5_i laws are properties of the solutions of a differential equation, symmetries are by definition 



properties of the differential equations and only in particular cases, namely in that of group- 
invariant solutions, also properties of the solutions of differential equations. Therefore, it is 
a standing question whether a discretization scheme that preserves numerically a property of 
a differential equation also improves the quality of the numerical solution of that discretized 
differential equation. 

The present paper is devoted to an investigation of this question and related problems ex- 
emplified with invariant discretization schemes for the linear heat equation. The heat equation 
is especially suited for this investigation as there already have been several studies devoted to 
the construction of invariant numerical schemes for this equation [1, 13, 34]. At the same time, 
in none of the existing works on that equations a deeper background analysis of the numerical 
properties (e.g. order of approximation or stability) of the developed schemes was given. 

There are several reasons why less attention has been paid so far on the numerical analysis 
of invariant schemes (with the exception of the work [10]). One of the reasons is that the field 
of invariant discretization schemes is still in its early stages, with new conceptual algorithms 



being developed only recently [2, 21, 22, 29]. Another reason is that invariant finite difference 
schemes generally require the usage of adaptive moving meshes, i.e. it is necessary to include a 
non-trivial mesh equation in the discretization problem. Moving meshes can lead to non-uniform 
grid point distributions and in the multi-dimensional case to non-orthogonal grids. The analysis 
of schemes on such meshes is considerably more difficult than that for related difference schemes 
on fixed, uniform and orthogonal meshes. Due to this second reason, most invariant numerical 
schemes so far have been constructed only for (l+l)-dimensional single evolution equations, 
as in the case of only one space dimension a moving mesh can be handled still with limited 
effort. Although we will be concerned with a (l+l)-dimensional equation in the present paper 
as well, the methods used subsequently can be employed in the multi-dimensional case without 
any substantial modification. 

The new approach we propose here is to use the invariant grid equation only for a single time 
step and then to interpolate the solution to the regular grid. The important observation is that 
this interpolation can be done in an invariant way, i.e. re-mapping the solution does not break 
the invariance of the scheme. At the same time, the possibility to re-map the solution of an 
invariant scheme to a regular grid is highly desirable as in the multi-dimensional case a freely 
evolving grid can cause severe numerical problems. Moreover, for realistic numerical models, 
as e.g. employed in weather and climate predictions, it will in general be impossible to use 
adaptive meshes as the discretization of the governing equations is only one part of such model. 
Other parts are related to the numerical data assimilation, i.e. the preparation of the initial 
conditions for the numerical model and this step usually involves the forecasting model itself. 
As the assimilation of the initial conditions cannot be done on an evolving mesh (because the 
data are given at fixed locations only) this at once renders invariant schemes on moving meshes 
highly impractical. Equally important, any realistic numerical model for a nonlinear system 
of partial differential equations has to contain subgrid-scale models, which mimic the effects 
of processes taking place at those scales that the numerical model is not capable of resolving 
explicitly [32, 33]. The construction of subgrid-scale models for non- resolved physical processes 
in general involves ad-hoc arguments and these arguments rely on the particular scale on which 
the unresolved processes take place. As a moving mesh locally changes the resolution and thus 
has an impact on what processes are explicitly resolved by a numerical model, subgrid-scale 
schemes have to be designed that can operate on grids with varying resolution. For realistic 
processes (which are usually not self-similar), this might be very hard to achieve in practice. 

All what was said above are strong objections against using invariant numerical schemes for 
multi-dimensional systems of differential equations on freely evolving meshes and thus, whether 
mathematically feasible or not, would only attract little practical interest. This is why other 
approaches should be sought that on the one hand allow one to retain the invariance group of 
a system of differential equations in the course of discretization and on the other hand yield 
schemes that are practical enough to avoid the above mentioned and related problems. The 
proposed evolution-re-mapping strategy we are going to introduce below may be considered as 
one such approach. 

The further organization of this article is as follows. The subsequent Section 2 includes a 
brief summary and some extensions on the various methods to construct invariant discretiza- 
tion schemes. In Section 3 the heat equation along with its maximal Lie invariance group G 
is presented. It is discussed which subgroup G 1 of G we aim to present when constructing in- 
variant numerical discretization schemes. The selection of G 1 is based on preserving the class 
of initial-boundary value problems we are focussing on. Section 4 contains the construction of 
an equivariant moving frame for G 1 along with a presentation of some lower order differential 
invariants of G 1 . In Section 5 a simple invariant discretization scheme for the heat equation in 
computational coordinates is found. The local discretization error of this scheme and related 
schemes is established in Section 6. In Section 7 we introduce the new idea of invariant interpola- 
tion schemes that can be used to map the numerical solution obtained from an invariant scheme 
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on a moving mesh back to a regular grid. The numerical analysis as well as some numerical 
tests for the schemes proposed in this paper are found in Section 8. The summary of this article 
is presented in the final Section 9. 



2 Construction of invariant discretization schemes 

The construction of invariant discretization schemes for differential equations can be seen as a 
part of the ongoing effort to turn the group analysis into an efficient tool for the analysis of 
difference equations, see e.g. the review article [24]. As of now, there are three main methods 
that are used to construct discretization schemes with certain invariance properties. 

2.1 Difference invariant method 

The first method was developed by Dorodnitsyn, see [1, 12, 13, 24, 34]. It operates with the 
infinitesimal generators v E Q of one-parameter symmetry transformations that span the maxi- 
mal Lie invariance algebra g of the system of differential equations under consideration. These 
generators are of the form 

v = ( j (x,u)d x j + r] a (x,u)d u <* = ((x,u)d x + rj(x,u)d u , 

where x = (x 1 ,...,x p ) and u = (u , ...,u q ) are the tuples of independent and dependent 
variables, respectively. Here and in the following, the summation convention over repeated 
indices is used. Rather than prolonging v to higher order derivatives of u with respect to x, 
which is standard in the symmetry analysis of differential equations [3, 25, 28], in this method 
the vector fields are prolonged to all the points of the discretization stencil, i.e. the collection 
of grid points which is necessary to approximate a given system of differential equation up to a 
desired order. This prolongation is of the form 



where X\ = (xj, . . . , x?) and U{ = (uj, . . . , uf), i.e. it is formally done by evaluating a vector field 
v at all the m stencil points zj = (xi,U{) and summing up the result. An example for such a 
prolongation is given in Remark 1 in Section 5. 

As a next step, the invariants of the group action are found by invoking the infinitesimal in- 
variance criterion [25], which in the present case is prv(J) = 0, whenever 1 = holds. Functions 
/ that fulfill this condition for all v E Q are termed difference invariants. 

Once all the difference invariants on the stencil space are found, one then attempts to assemble 
these invariants together to a finite difference approximation of the given system of differential 
equations. By construction, this procedure guarantees that the resulting numerical scheme is 
invariant under a symmetry group that is isomorphic to the symmetry group of the original 
system of differential equations. 

The main drawback of this method is that it might be hard to find a good set of difference 
invariants that allows one to approximate differential equations in the multi-dimensional case. 
The problem is, as discussed in the introduction, that invariant schemes generally require the 
usage of moving and/or non-orthogonal grids. Formulating grid equations using difference in- 
variants and approximating finite difference derivatives in an invariant way on the resulting grids 
can be rather challenging in higher dimensions and thus limited the application of this method 
to the case of single (1 + l)-dimensional evolution equations, although this limitation is merely 
of technical than conceptual nature. 
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2.2 Invariant moving mesh method 

The invariance of finite difference schemes under the maximal Lie invariance groups of several 
physically relevant time-dependent differential equations, requires the usage of moving meshes. 
This is true both for the finite difference method discussed in the previous Section 2.1 and the 
moving frame method to be discussed in the next Section 2.3. This kind of mesh adaption in 
which the number of grid points remains constant throughout the integration is referred to as 
r-adaptivity in the field of adaptive numerical schemes [7, 20]. 

The standard strategy to handle r-adaptive meshes is to regard the grid adaption as a time- 
dependent mapping from a fixed reference space of computational coordinates to the physical 
space of the independent variables of the differential equation, i.e. x = for £ = (£ , . . . ,£ p ) 
being the computational variables. Without loss of generality, we assume that = r = t is the 
time variable. Expression of the dependent variables u in the computational space is achieved 
by setting k(£) = u(x(£)). 

The significance of the computational coordinates is to provide a reference frame that remains 
fixed and orthogonal even in the presence of grid adaption in the physical space of coordinates. 
In the course of discretization, the variable £ labels the position of the grid points in the mesh 
and this labeling stays unchanged during the mesh adaption. Thus, the computational variables 
can be interpreted physically as Lagrangian coordinates and their invariance under the motion 
of the grid is equivalent to the identity of fluid particles in ideal hydrodynamics. 

Because by construction the grid remains orthogonal in the ^-coordinates, the usual finite 
difference approximations for derivatives can be used in the space of computational variables. 
This simplifies both the practical implementation of the discretization method and the numerical 
analysis of the resulting schemes. 

The expression of the initial physical system of differential equations in terms of computa- 
tional variables leads to a system of equations that explicitly includes the mesh velocity x T , 
which is yet to be determined in order to close the resulting numerical scheme. The strategy for 
determining the location of the grid points at the subsequent time level follows from the equidis- 
tribution principle, which in its differential form is {px^)^ = 0, where p is a monitor function 
that determines the areas of grid convergence and divergence. The drawback of this approach 
is that the equidistribution principle completely determines a grid only in the case of one-space 
dimension. For higher-dimensional problems, equidistribution has to be combined with heuristic 
arguments, see [20] for more details. 

The invariance of the initial differential equations is brought into the scheme by adequately 
specifying the monitor function p. In [6] it was pointed out that using monitor functions that 
preserving the scale-invariance of a differential equation is especially relevant in cases where 
the equation is capable of developing a blow-up solution in finite time, see also [7, 8, 19]. This 
finding is generalized upon requiring that the monitor function is chosen in a manner such 
that the equidistribution principle is invariant under the same symmetry group as the original 
differential equation. For a number of symmetry groups this appears to be possible, see [2] for 
an example. 

The invariant moving mesh method was recently extended in [2] . The idea of this extension is 
to transform the initial system of differential equations to the space of computational coordinates 
and to determine the form of the action of the symmetry transformations in the computational 
space. The equations in the computational space are then discretized in such a manner that the 
resulting scheme transforms according to the discrete version of the transformation laws found. 
The main advantage of this approach is that it allows one to retain an initial conserved form 
of the system of differential equations and thus to numerically preserve certain conservation 
laws in the invariant schemes. This is relevant as preserving conservation laws in the course of 
invariant numerical modeling is yet a pristine problem. An exception to this is the discretization 
of equations that follow from variational principles, which, if done in a proper way, can lead 
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to the simultaneous preservation of both symmetries and associated conservation laws, owing 
to the discrete Noether theorem. See, e.g. [5] for an example of such an invariant Lagrangian 
discretization. 

2.3 Moving frame method 

The third method is the most recent one [10, 11, 21, 22, 29]. It relies on the notion of equivariant 
moving frames and their important property to provide a mapping that allows one to associated 
to a given function an invariant function in a canonical way. As we will mostly work with this 
method in the present paper, we describe it in greater detail here. We collect some important 
notions on moving frames below, a more comprehensive exposition can be found in the original 
references [9, 14, 15, 26, 27, 29]. 

Definition 1. Let G be a finite dimensional Lie group acting on a manifold M. A (right) 
moving frame p is a smooth map p: M — >■ G satisfying the equivariance property 

pig ■ z) = p{z)g~ 1 , (1) 

for all z G M and g G G. 

Theorem 1. A moving frame exists in the neighborhood of a point z G M if and only if the 
group G acts freely and regularly near z. 

Freeness of a group action means that z = g ■ z = z for all z G M only holds for g being the 
identity transformation, which implies that all the group orbits have the same dimension. Here 
and throughout the paper, a tilde over a variable denotes the corresponding transformed form 
of that variable. Regularity of a group action requires that there exists a neighborhood for each 
point z G M, which is intersected by the orbits of G into a pathwise connected subset. 

When a group G does not act free on M, it can be made free upon extending the action 
of G to a suitably high-order jet space J n = J n (M,p) of M, < n < oo. Locally, the nth 
order jet space of a p-dimensional submanifold of M has coordinates z^ = (x,u^), where 
as in the previous subsections x = (x ,...,x p ) are considered as the independent variables, 
u = (it 1 , . . . , u q ), q = dim M — p, are the dependent variables and collects all the derivatives 
of u with respect to x of order not greater than n including u as the zeroth order derivatives. In 
practice, the prolongation of the group action of G on J n is implemented using the chain rule. 

Moving frames are determined using a normalization procedure. The steps to find a moving 
frame for a group action G are the following: (i) Define a cross-section to the group orbits. A 
cross-section C is any submanifold C C M of complementary dimension to the dimension r of 
the group orbits, i.e. dimC = dim M — r, that intersects the group orbits once and transversally. 
Usually coordinate cross-sections are chosen in which some of the coordinates of M (or of J n 
if the group action is not free on M) are set to constants, i.e. z 1 = c , . . . , z r = c r . (ii) The 
algebraic system z 1 = g ■ z 1 = c , . . . , W = g ■ z' r = c r is solved for the group parameters 
g = (ei, . . . , e r ). The resulting expression g = p{z) is the moving frame. 

Moving frames can be used to map any given function to an invariant function by a procedure 
called invariantization. 

Definition 2. The invariantization of a real- valued function / : M — )■ R using the (right) moving 
frame p is the function which is defined as i{f){z) = f(g ■ z)\ g=p ( z \ = f(p(z) • z). 

That the function constructed in this way is indeed invariant follows from the equivari- 
ance property (1) of the moving frame p, 

<f)(9 ■ *) = f(p(9 ■ z)g ■ z) = f(p(z)g- 1 g • z) = f(p(z) ■ z) = t(f)(z), 
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which boils down to the definition of an invariant function /, i.e. I(g ■ z) = I(z). In practice, a 
function f(z) is invariantized by first transforming its argument using the transformations from 
G and then by substituting the computed moving frame for the group parameters. By definition, 
an invariant that is defined on the jet space J n is called a differential invariant. 

Moving frames can also be constructed on a discrete space. In a finite difference approxima- 
tion the coordinates on J n , i.e. derivatives, are approximated using a finite set of points, and all 
the points needed to approximate the derivatives arising in a system of differential equations are 
the points of the stencil introduced in Section 2.1. Because most of the interesting symmetries 
of differential equations that are broken in standard numerical schemes require the usage of 
non-orthogonal spatial-temporal discretization meshes, it is beneficial to both regard x and u as 
the dependent variables and the computational variable £ as the independent variables as was 
discussed in the previous Section 2.2. 

Sampling the tuples from the extended computational space = {(f,,z)} at discrete points, 
i.e. at = zf), one can introduce the space M| n = {(w\, . . . , w n ) | & 7^ £j for all i / 

j}, where Wi = which can be identified as the joint product space of stencil variables. 

Because the identifier £j of points z are required to be unique, each element of M^ n only includes 
distinct grid points in the physical space of equation variables. The dimension of the space M^ n 
depends on the number of independent and dependent variables in the system of differential 
equation and the desired order of accuracy of the approximated derivatives. 

It is possible to carry out the construction of the moving frame on M£ n , i.e. to define the 
moving frame by an equivariant mapping p^ n : M^ n — > G, where G acts on M^ n by the product 
action, g ■ (w\, . . . ,w n ) = (g • w%, . . . ,g • w n ). Note that the extension of the group action to 
the computational variables £ is trivial, i.e. £ = g ■ £ = £, they remain unaffected by G, see [2]. 
The compatibility between the moving frame p^ n and the moving frame p on the space M (or 
an appropriate jet space J n ), i.e. that p^ n — > p in the continuous limit is assured provided 
that the cross-section defining the moving frame p^ n in the continuous limit converges to the 
cross-section defining the moving frame p. Once the moving frame is constructed on the discrete 
space M^ n of stencil variables, it can be used to invariantize any numerical scheme expressed in 
computational coordinates. This will be explicitly shown in Sections 5 and 6, where we construct 
invariant schemes for the heat equation. 

It is essential that the construction of the moving frame on the grid point space is carried 
out in terms of computational rather than in physical coordinates. This can be illustrated by 
the following simple example. 

Example 1. The Laplace equation u xx +u yy = is, inter alia, invariant under the one-parameter 
group of rotations SO(2), x = x cos e—y sine, y = x sine+y cose. Let us obtain the moving frame 
p for this group action from the normalization condition u x = 0, i.e. we determine the moving 
frame on the first jet space J 1 (M, 2), p = p(x,y,u,u x ,u y ). Prolonging the transformations 
from 50(2) to the derivative u x leads to u x = u x cose — u y sine and thus the moving frame is 
e = &rct&n(u x / u y ) . 

Let us now attempt to find the moving frame from the discrete normalization condition 
u x = 0, where we choose centered differences to evaluate the discrete derivative. According 
to the theory, the resulting moving frame is compatible with the moving frame p as in the 
continuous limit u x — > u x and thus the cross-sections are compatible. However, evaluating u x 
in the nai've way = {ui+\ — — Xj_i), we fail as 

„d = g*±i ~ = o 

x (x i+ i -Xi-i) cose- (y i+ i - yi-i) sine 

cannot be solved for the group parameter e. On the other hand, setting u = u{x(£}, £ 2 ), £ 2 )) 
and expressing u x in terms of the computational variables £ 2 , the normalization u x = reads 
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as 



u d 



uji y^ 2 - u^ 2 y^ u|i (xj 2 sin e + cos e) - (x^ sin e + cos e) 



/yd /i ,d rfd n id /yd /i ,d rv*d n ,(1 

This expression can be solved for e and it gives 
e = arctan — ^ — ^ ^ — ^- = arctan — 7 



which in the continuous limit goes to e = &ictan{u x / u y ) as required. 



3 Lie symmetries of the heat equation 

The one-dimensional linear heat transport equation is 

u t - u xx = 0, (2) 

where we scaled the thermal diffusivity v to 1, i.e. Eq. (2) is in non-dimensional form. 

The heat equation (2) admits the following infinitesimal generators of one-parameter groups, 
which generate the maximal Lie invariance algebra g of Eq. (2): 

d t , d x , ud u , 2td t + xd x , 2td x - xud u , 
At 2 d t + Atxd x - (x 2 + 2t)ud u , a(t,x)d u , 

where a is a solution of Eq. (2), see e.g. [25]. These vector fields generate (i) time-translations, 
(ii) space translations, (iii) scalings in u, (iv) scalings in t and x, (v) Galilean boosts, (vi) 
inversion and (vii) the superposition principle symmetry. 

In what follows we will focus our attention on numerically preserving the symmetries gen- 
erated by the subalgebra rj 1 spanned be the first five operators. We do not require the sixth 
operator to be preserved as this operator does not map a given initial-boundary value problem 
of Eq. (2) to another initial-boundary value problem of Eq. (2), i.e. it does not preserve the form 
of a given class of such problems for the heat equation. In particular, the inversion symmetry 
does not send an initial-value problem to another initial-value problem. See [2] for a further 
discussion about this problem. We also do not require to preserve the linearity operator here 
by construction. At the same time, as will be shown in Section (8), the numerical schemes we 
propose in this paper preserve the linearity property up to the discretization error expected. 



4 Moving frame and differential invariants for the heat equation 

We determine the moving frame for the subgroup G 1 of transformations associated with the 
subalgebra g . Transformations of G 1 are of the form 

t = e 2£4 {t + x = e £4 (x + £ 2 + 2£ 5 t), u = e E ^ x ~^u. (4) 

Because the group action of G l is not free on M = {(t, x, u)} we have to carry out the moving 
frame construction on a suitable high-order jet space. In the present case, the group action of 
G l becomes free on J 1 = J 1 (M, 2). Thus, it is necessary to extend the transformations (4) to 
derivatives of u with respect to t and x. 

Using the chain rule we can compute the transformed operators of total differentiation, which 
read as 

B i = e- 2£4 (D t -2e 5 B x ), B x = e~ £4 B x , 
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where D x = d x + u x d u + u tx du t + u xx d Ux + ... and D t = d t + u t d u + u tt d Ut + u tx d Ux + ... denote 
the usual operators of total differentiation. With the transformed total differentiation operators 
at hand, it is possible to compute the transformed partial derivatives of u with respect to t and 
x. The transformation rules for the lowest order derivatives are 



_ 2£4+£3 _ e5X _ e 2^ _ 2£s ^ + £ 2 u)) = e - £i+£3 - £5X - e lt^ x _ 



-5 

=■2 



u xx = e- 2£ ^- e ^\u xx - 2e 5 u x + s 2 5 u). 

In fact, for the construction of the moving frame already the knowledge of the first order deriva- 
tives is sufficient. 

We determine the moving frame for the five-parameter group of transformations of the 
form (4) by using the following normalization conditions which determine a valid cross-section 
to the group orbits of the prolonged action of G 1 on J 1 , 

t = 0, x = 0, u = 1, ut = 1, u x = 0. (5) 

The moving frame is computed by taking the transformed form of the normalization conditions, 
t = 0, x = 0, u = 1, tt 4 - = 1 and u x = and by solving the resulting algebraic system for the group 
parameters £i, . . . , £5. The result of this computation is the following moving frame g = p(z^), 



2 

u x \ I u x u x 



£\ = — t, £2 = — [x , £3 = — I hxu — x t 



u / \ u u 2 



u t u 2 x u x 



(6) 



£4= In W 1, £5 



u u 2 u 



With the moving frame at hand, we can invariantize any of the partial derivatives of u with 
respect to t and x and thus obtain a complete set of differential invariants for the subgroup G 1 
of the maximal Lie invariance group of the heat equation. As an example, invariantizing the 
derivative u xx , i.e. computing i(u xx ) = (g ■ u xx )\ g=p ^ z (i)^ we produce the differential invariant 



2 2 
j _ uu xx — u x 

UUf — u 2 



Setting / = 1, we obtain 

u(u t - u xx 



UUt — u 2 



0. 



which yields the original heat equation expressed in terms of differential invariants. This re- 
expression of a differential equation using differential invariants is know as the replacement 
theorem [9]. 



5 Invariant discretization of the heat equation 

The discretization of Eq. (2) cannot be done on a fixed, uniform grid. To see this, let us check 
the transformation behavior of the grid equation — xf = 0, which is the definition of a 
stationary grid, under the transformations (4). Here and in what follows, a hat stands for a 
variable evaluated at the time t n+1 = t n + At. This yields 

_„ +1 _ -„ = e e^ x n+l _ x n + 2e b (f+ 1 -f 1 )), 

which is only zero in the case when £5 = 0. Stated in another way, a discretization on a fixed 
grid can at most preserve the symmetry subgroup of G, which is generated by the first four 
elements of the maximal Lie invariance algebra q of the heat equation (3). 
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Thus, the discretization of (2) preserving G 1 will require the usage of non-stationary grids. 
For this reason it is convenient to express (2) in terms of computational coordinates initially, 
i.e. we set ii(r, £) = u(t, x(t, £)), where £ is the single spatial computational variable and r = t. 
The heat equation in this set of coordinates reads 



x- 



3 



Uc 



0. 



(7) 



So as to find the invariant discretization of the heat equation in the form (7), we determine the 
moving frame in the space of stencil variables M^ 4 using the discrete analogs of the normalization 
conditions (5), but expressed in terms of computational coordinates. 



For the sake of convenience, we introduce the notation, h? 



At 



r n+l 



xf, fi 



r n . The discretization stencil we use is depicted in Fig. 1. 



■ + At 



(a 



Figure 1. Stencil for the invariant discretization scheme of the heat equation. 
The appropriate normalization conditions for a compatible moving frame p^ 4 are 
r n = 0, x™ = 0, < = 1, 4 = 1, u d x = 0, 

where 



(8) 



.71+1 



d U i+1 



U 



i-1 



a 



At 



h+ + fi- 



ll', 



t+l 



H-l 



h+ + fl- 



are the discretizations of the time and space derivatives expressed in computational coordinates 
and x d = xf)/ At is the discrete grid velocity. Replacing the single equations in the above 

normalization conditions by their respective transformed expressions and solving the resulting 
algebraic system for the group parameters we obtain the following moving frame on the space 
of stencil variables M^ 4 , 



ei = -T n , e 2 = -(x?-2T n (lnu x ) 



£3 



e 4 



1, /exp {-AT{x d T {\nu x ) d + ((lnu^) 2 )) 
in 1 



In up 



n+l 



x?(lnu x ) d -T n ((lnu, 



\d\2 



u?At 



e 5 = (lnit x ) 



d 

X) 1 



(9) 



where we introduced (\nu x ) = (lnu^ +1 — lnu"_ 1 )/(/i + + h ). This moving frame is compatible 
with the moving frame (6) in that it converges to (6) in the continuous limit A£ — > and 
At -> upon using h + = x f A£ + x#{A£) 2 /2 + 0(A£ 3 ), h~ = x ? A£ - x f€ (A£) 2 /2 + <3(A£ 3 ), 



-71+1 



x? + x t At + 0(At 2 ), uf +1 = uf + u^At + u a (A£) 2 /2 + 0(A^ 



i, <*i_i 



u K (A0 2 /2 " 0(A^ 3 ), and u™ +1 = < + ii T At + 0(Ar 2 ). 

The moving frame (9) can now be used to invariantize any non-invariant finite difference 
discretization of (7) on M^ 4 . To illustrate this, we start with the standard FTCD scheme 

1 ' ' " ' (h+-h-)4)=o, 



(h+ + h-) 2 



i+l 



+ u, 



i-1 



2u? 
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and apply the invariantization map to it. This is done by first replacing all the single expressions 
by their respective transformed expression and substituting the moving frame (9) for the arising 
group parameters. The result of this procedure is the following invariant scheme 



_ exp (-At {x d T (lnu) d x + ((ln^) 2 )) < +1 - < 

At 

h+/(h++h-) , n \ h~ /(h + +h~) 



u i-lJ \ u i-l 



(h+ + h-) 2 



(10) 

0, 



Again, it can be checked that the above scheme (10) indeed converges to (7) in the limit of 
A£ —7- and At — > 0. This will be shown explicitly in Section 6. 

So as to complete the scheme (10) it is necessary to derive and equation for which 
is the only ingredient missing in (10). There are different ways to determine a grid equation, 
such as using the equidistribution principle as outlined in Section 2.2. The problem with this 
strategy in the present case is that while it might be beneficial from the numerical point of view, 
it might not be easy to obtain an invariant discretization of this principle which does not lead 
to a fully coupled equation-grid system, i.e. it can happen that the grid equation includes both 
u™ and u™ +1 . While this coupling is not a problem in the one-dimensional case, it can lead to a 
severe restriction of the applicability in the multi-dimensional then solving the coupled 

equation-grid system might be too expensive and other strategies must be pursued. 

We circumvent the aforementioned coupling problem by using a grid equation that does not 
involve u™ +1 , such as 

which can be checked to be invariant under the subgroup G 1 of the maximal Lie invariance 
group of the heat equation, see Remark 1. This is the same grid equation as was used in [13]. 
In the continuous limit, this grid equation goes to 

2 

x T = (lnu)c. 

H 

Remark 1. While the invariantization algorithm guarantees that the scheme (10) is indeed 
invariant under the subgroup G 1 of the maximal Lie invariance group of the heat equation, this 
invariance can be checked in a straightforward fashion using the infinitesimal invariance criterion 
as invoked in the Dorodnitsyn method discussed in Section 2.1. Let us recall that this criterion 
states that an invariant I of a group action is annihilated by the associated infinitesimal gener- 
ators, i.e. v(I) = when 1 = holds for all v G fj. Because in the present case, the invariants 
are defined on the stencil space with coordinates T n , At, xf, x™_ l5 uf, it£_ ij , 
we have to prolong the operators of q accordingly. The prolongations of the first five operators 
of (3) to the variables of the stencil are 

d T n, 8 x n + 8 x n + d x « + 8 n + 1 , U™ d u ™ + uf, 1 8 U " +ti™i9 n ™ + U" + 1 8 n + 1, 

2r n d T n + 2Ard Ar + x?d x n + x? +1 8 x?+i + x?_ 1 9 a . ? _ 1 + x? +1 d x n+i , 2r n (d x? + d x?+i + 
d xf J + 2(r n + Ar)d x7+l - xfu?d K - " " x 7 +l < +1 V 1 ' 

see e.g. [12] for more details. It can be checked that pr v(5) = and pr v(M) = hold on S = 
and M = for all the prolonged infinitesimal generators and thus S = is a proper invariant 
numerical scheme and M = an invariant grid equation. 
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6 Numerical properties of the invariant scheme 



In this section we investigate the numerical properties of the scheme (10) and related schemes. 
We start our consideration with the estimation of the local truncation error of the scheme. The 
study of this question is relevant because so far little is known about the relation between the 
order of a non-invariant scheme and its corresponding invariantized scheme. 

The discretization of the heat equation in computational coordinates (7) can be formally 
represented as 

u "- x "|-(^( 4 " ^N) =0 ' (12) 

where in the present case we assume that derivatives are approximated with the aid of a standard 
FTCD scheme. More general schemes will be considered after the order of the invariant FTCD 
scheme is established. 

Theorem 2. The order of the invariant scheme (10) is the same as the order of the scheme (12), 
namely first order in time and second order in space, provided that an Euler forward step and 
second order centered differences are used to approximate the time and space derivatives arising 
in both the differential equation (12) and the normalization conditions (8), respectively. 

Proof. Invariantizing the scheme (12) using the normalization conditions (8) leads to 

where i(f(z)) denotes the invariantization of the function f{z). By definition, invariantization of 
a function f(z) means to transform the argument z and plug in the moving frame for the group 
parameters. In the present case, the transformed form of the last expression can be written as 



4 e £3-£5Z-£5*-2£4 
r d\2 



1 " ~ , it (e~ £5h X+i + e^ h ~uU ~ 2<) = 0. 



Using the normalization condition m" = 1 we obtain that 

uf = 1 = e e3-e 5 x-e?J u n 

and thus the last expression can be recast as 

e2£% " " {h+ + h - ) 2 ^ £5h+ ^ + e£5 ""<-i " 2u i) = 0- (14) 

Let us now determine the local discretization error in the parameter e§. The respective moving 
frame component is 

_ ln< +1 -ln<„ 1 
£5 ~ h+ + h~ ' 

which, upon using h+ = x ? A£ + x^(A^) 2 /2 + 0(A£ 3 ), h~ = x 6 A£ - x^(A£) 2 /2 + 0(A£ 3 ), 
<+i = < + u 5 A£ + u €f (A£) 2 /2 + 0(A£ 3 ), u ^ = u ^ - Ui :A^ + u^(A0 2 /2-O(A^) expands to 

e 5 = -^ + 0(Af). (15) 
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Substituting £5 into the second term of Eq. (14) and expanding the exponential functions in the 
same term into Taylor series, we obtain after some rearranging 

e 2£4 < - x 2 Ae l 0{Ae) + <-i - K - £ 5 (x 5 Ae« +1 - ujLi) + 

l -x^Af(v^ +1 + + ^^|A^ 2 (nr +1 + 11SL1) + 0(A^ 4 )) = 0, 

which can be further simplified to 

e ^ u n _ J_ ( _ !| _ m j + 0(A ^ 2) = (16) 

It now remains to expand the first term in Eq. (16). The moving frame component for £4 in (9) 
can be recast as 

' — x 1 } lnu^ — ln*u™_i /lnii™i_i — hiit? 



e^pj ^ At -- ^ + ( "7- ' 1 1 < +1 

At 




Using = xf + rr r At + 0(Ar 2 ) and u™ +1 = uf + u t At + 0(Ar 2 ) and again expanding the 
exponential function into a Taylor series, we derive 



e 2 ^ = Ut - Xt ^L-W + (At, 
X£ x^ 



Plugging this into Eq. (16) we arrive at 

u T - x^ - \ (u^ - —uA + 0(At, Af) = 0, 

X^ X p \ X^ J 



which completes the proof of the theorem. □ 

A more general statement is the following one: 

Theorem 3. The order of the spatial discretization of an invariant finite difference scheme of the 
heat equation in computational variables preserves the order p G N of the spatial discretization of 
the associated non-invariant finite difference scheme provided that centered differences of order 
p are used to approximate both the derivatives in the heat equation and in the normalization 
conditions (8). 

Proof. In view of the general form (13) of the invariantization of the scheme (12), we have to 
study the effect of invariantization on the discretizations of the terms x^ and u^. 

The invariantization of (xf) 2 = (x^ + 0(A£ P )) 2 is simply t((x|) 2 ) = e 2£4 (xj + 0(A£ P )) and is 
of the same order p if, as required, the moving frame component £4 stems from the approximation 
of u~: = 1 using pth order accuracy and thus only includes approximations of derivatives with 
that accuracy. 

Let us now investigate the invariantization of discretizations of u«. The general form of a 
centered difference approximation of of even order p is 

1 P/2 
i=-p/2 
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with coefficients c 2 pd = 2^/j, jeA = {-p/2, ...,-1,1,... ,p/2}, c 2 pfi = -2 £f=i V* 2 and 

x {-iy+\ P /2r 

™ j(p/2 + j)\(p/2-j)V 
for j £ A and c^q = are the coefficients from the pth order approximation of u%, i.e. 

1 P/2 

see [16] for a discussion of the algorithm for finding the weights ct • in higher-order finite differ- 
ence approximations of the feth derivative of u. The invariantization of is 

1 P/2 

i=-p/2 

or, upon using the normalization condition -u™ = 1, 
1 1 P/2 

i ( n «) = 7 x 72^ E c^.exp(-e 5 Ax j )^, (17) 
* 4 i=-P/a 

where we expand 

k=l s 1=0 s 

Using the expansions for Axj and u™, the expression (17) can be expanded and rearranged in 
powers of jA£ in the form 

1 1 P/2 DO 

4(u « ) = ap^ E £(-i)^o-A6 fc , 

^ 4 j=-p/2fc=0 



where 

1 

^2 = - (it# - £5(2x 5 u € + x^Ui) + elxju^) , 



and the expressions for Aj., k ^ 2 are of no importance. The proof is completed upon substituting 
for £5 the corresponding moving frame component (which is of order p when the normalization 
u^. = is approximated with pth order accuracy) and by noting that 



P/2 



E ^ 



for ke {0,1,3,... p + 2, 2n}, re e N 

2 for k = 2 

Cfc 7^ else 



j=-p/2 

where the precise values of the constants c k follow from evaluating the respective sums. □ 

The scheme (10) is only of first order in time r = t. To construct a scheme that is second 
order in time, we can start with a non-invariant scheme (12) and discretize the time derivative 
u% with second order accuracy, e.g. we set uf = (u™ +1 — u™~ )/ (2 At), where re" -1 is the value of 
u at the previous time step r n_1 = r n — At. It is now necessary to check whether invariantizing 
this leapfrog discretization leads to an invariant scheme that is also second order in time. 
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Theorem 4. Invariantization of the scheme (12) in which a leapfrog step and second order 
centered differences are used to approximate the time and space derivatives, respectively, leads to 
an invariant scheme that is both second order in time and space provided that the normalization 
conditions (8) are approximated also using discretizations that are of second order. 

Proof. To prove this theorem, it is only necessary to establish the order of the first term in 
Eq. (16). We proceed in an analog manner as in the proof of Theorem 2, i.e. we discretize the 
normalization condition u~. = 1 but now with second order accuracy. This yields 

, _ e * 5 / ega^AT- e2AT„n+l _ „e 5 x T Ar+e 2 5 A-rn-1 



~a = _ -^ar-qaT nti _ n-i = j 

2e 2£ 4Ar V 4 1 J 

Using the normalization condition uf = 1 as before and expanding the exponential functions we 
derive 

e 2£4 < = ^ « +1 " <~ X " Ar(e 5 x T + ef)^ 1 + u?- 1 )) + 0(Ar 2 ) 
and upon noting that u™ +1 + = 2u™ + 0(Ar 2 ) we obtain 



Xg uf 



where we have substituted the expression (15) for £5. Plugging this result into Eq. (16) completes 
the proof of the theorem. □ 

The actual form of the resulting invariant leapfrog scheme is the following 

exp (-At {x d T (lnu) d x + ((ln^) 2 )) < +1 - exp (At (^(ln^ + ((ln^) 2 )) u?" 1 



S 



2Ar 

(h+ + h-) 2 ' 

(18) 

where xf = (x" +1 - xf )/At and = (xf - x" _1 )/At. 

Higher order in time schemes can be constructed upon invariantizing multi-stage schemes. 
Combining this result with the result established in Theorem 3 we have found the following: 

Corollary 1. Invariantizing a non-invariant finite difference scheme for the heat equation in 
computational coordinates preserves the spatial and temporal order of the initial non-invariant 
finite difference scheme provided that the normalization conditions for the moving frame are 
discretized with the same order as the respective derivatives in the non-invariant finite difference 
scheme and centered differences are used. 



7 Invariant interpolation schemes 

A common property of invariant numerical schemes for evolution equations possessing a nontriv- 
ial maximal Lie invariance group is that it is not possible to use a fixed, orthogonal discretization 
mesh. The continuous evolution of the mesh, if not handled properly, can lead to several unde- 
sirable properties, such as an overly strong concentration of grid points in certain regions and 
therefore too poor a resolution in other parts of the integration domain. The problem gets worse 
in the multi-dimensional case, where mesh tangling or strongly skewed meshes can occur. But 
even if the mesh movement can be managed in an optimal way there are various problems for 
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which it is impossible to use continuously adaptive grids. An example for this are practically all 
models that are in operational use in weather and climate prediction. These models employ so- 
phisticated data assimilation strategies and are coupled to subgrid-scale parameterizations that 
aim to mimic the effects of unresolved processes on the grid scale variables. Attempting to make 
use of data assimilation or parameterization schemes on moving meshes is not only a technical 
problem that would cause a significant computational overhead compared to standard schemes 
but also a conceptual challenge for it is unclear on how to design parameterization schemes that 
can operate on grids with varying resolution. In order to promote the ideas of invariant nu- 
merical discretization schemes beyond their application to simple evolution equations it is thus 
instructive to study possible ways of overcoming the limitations imposed by the requirement of 
using moving meshes. 

One straightforward idea is to use invariant schemes on fixed (i.e. non- invariant grids). As 
was shown in [29] this can lead to improved numerical solutions compared to non-invariant 
integrators, while still being excelled by the results that can be obtained if using completely 
invariant schemes. On the other hand, if moving (invariant) meshes are not tractable for a 
particular class of problems, preserving the invariance of a system of differential equations at 
least for the discretization of this system itself might be a possible trade-off to take. 

Another idea is to use an evolution-re-mapping strategy, which will be proposed in the fol- 
lowing. This concept relies on using the invariant scheme for the system of differential equations 
together with the invariant mesh equations but the numerical solution is mapped back to a reg- 
ular mesh after each time step. A similar strategy has proven successful in the semi-Lagrangian 
time integration [31]. 

In the present case, the re-mapping step can be practically realized by using interpolation. 
Obviously, any standard interpolation scheme can be used to map the numerical solution 
defined at x™ +1 to the uniformly spaced £-grid. This, however, most likely breaks the invariance 
of the numerical scheme as a whole and so the question arises whether it is possible to accomplish 
the interpolation step in a symmetry-preserving fashion. 

In the following we discuss two possible ways of formulating invariant interpolation schemes, 
both of which can be used for finding interpolations that allow re-mapping of the numerical 
solution of a system of differential equations on a moving mesh to a Cartesian, equally-spaced 
grid. These ways are the invariantization of non-invariant interpolation schemes with moving 
frames and the construction of interpolations using difference invariants. 

Invariantization of interpolation schemes. The moving frame constructed in the course of 
invariantizing a finite difference scheme can also be used to invariantize a certain interpolation 
method. We exemplify this idea by invariantizing the formula for linear interpolation, 

where y G [s™ +1 , x i"+x ]■ The invariantization of this expression using a moving frame associated 
with G 1 yields the invariant interpolation formula 

< +1 (j/) = U? +1 + (y - < +1 ) Sj| - U ill , ^ +1 =exp((lnn £ ) rf (y-^ +1 ))< +1 . (19) 

x i+l x i 

Note that we have used a slightly different moving frame for the invarianization as we have used 
for invariantizing the finite difference scheme for the heat equation. Specifically, this moving 
frame is constructed by replacing the normalization condition u% = with ui = 0. The reason 
for this is that the moving frame used earlier yielded £5 = (lnu x ) d , i.e. it involves the solutions 
of Uj at the time step r n rather than at r n + At. Irrespectively of what normalization is 
used, both interpolations are invariant. Setting y = £j in the above interpolation formula yields 
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on the regular computational grid. Note that the interpolation (19) is consistent in that 
= < +1 and < +1 (x^) = lift 1 . 
In a similar manner more sophisticated interpolation schemes can be invariantized as well. 
In the following, we will use the invariantization of the quadratic interpolation formula. Usual 
quadratic interpolation is based on the expression 

n"+ 1 ( y )=<+%_ 1 ( y )+<+ 1 L J (y)+<+ 1 L m ( y ), L 3 (y) = \{ ^~^_ k n+l , (20) 

fc=»-l X j X k 

where y 6 Wi—i j^i+il an d Lj(y) are the Lagrangian interpolation polynomials. Invariantizing 
this formula using the same moving frame as above we get 

u n+l {y) = U^L^y) + U^Hy) + U?+ l L l+1 (y), (21) 

where Ui = exp((lnuf) d (x — x™ +1 ))u™ +1 , as in the case of the invariant linear interpolation (19). 
Numerical examples using the invariant quadratic interpolation will be given in Section 8. 

Interpolation using difference invariants. The moving frame on the grid point space allows 
invariantizing the elementary variables x™ and uf, which yields the system of joint invariants. 
In the continuous limit these invariants take the normalization values chosen for x and u to 
construct the usual moving frame p [26]. On the other hand, on the discrete space M^ n we only 
normalize one xf, i,n fixed, among all the grid points xf and the analog statement is true for 
the associated values uf. This means that the joint invariants i(xf) and i(uf), I ^ i, k ^ n, are 
nontrivial and can be used to assemble invariant interpolation schemes. 

In the present case, while we have normalized uf = 1 in the course of constructing the moving 
frame p* 4 , we are free to use the moving frame to invariantized any uf where I ^ i,k ^ n and 
this will yield a proper (nontrivial) invariant on the discrete space M^ 4 . As above, we again 
recompute the moving frame for G 1 by replacing the normalization conditions uf = 1 and 
u% = with = and ui = 0, respectively, which yields new expressions for the moving 

frame components of £3 and £5 given by 

£ 3 = -((ln< +1 - x? +1 (lnuz) d - r n+1 ((lnu £ ) d )) 2 ), £5 = (ln^) d . 

Using this modified moving frame, we then invarianize the variable u n+1 (£j), which is the sought 
value of u at the point (r n+ , £j) of the computational domain. This invariantization yields 

^ n+1 (&)) = ^Pexp((lnn £ ) d (x™ +1 -&)). 

Because in the continuous limit the invariantization of u n+l (£j) must reproduce the normalization 
condition u = 1, we restrict the difference invariant to the manifold t(it n+1 (£i)) = 1. The 
invariant interpolation is thus 

n" +1 fe) = exp((lnn £ ) d (£ 4 - ^ +1 )K +1 (22) 

and it is again consistent as u n+1 (x^ +1 ) = u™ +1 . More accurate interpolations could be con- 
structed by combining the invariants t>(xf) and t(uf) in a suitable way. 

The advantage of the interpolation methods introduced in this section is that they are invariant 
under the group G 1 , i.e. using these interpolation formulas to map u™ +1 back to the £-grid does 
not break the invariance of the numerical schemes for the heat equations, while still allowing 
one to use a regular grid and thus avoiding the complications that moving meshes impose on 
the applicability of symmetry-preserving finite difference discretizations. 
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8 Numerical verification 



In order to verify the accuracy predicted above for the various schemes proposed, we set-up the 
following problem. On a periodic domain x £ [0, 2n], consider 

Ut = Uxxi 

u(x,t = 0) = sin (2; - 1) + 2. 

On a sequence of grids with N £ {2, 4, 8, • • • , 256}, the number of grid points, we compute the 
error in the maximum norm between the numerical solution and the exact solution at t = 1. 
The time step dt is taken as proportional to h 2 , h = h + = hr , in all simulations. 

In each of the following figures, we plot the reference line corresponding to O (/i 2 ) in black, 
and the error as a red dotted line where, 

||.E|| L = max \u (x, 1) - Uexact (x, 1)1 . 
x&[0,2w] 

Note that for all approaches described below we expect a second order convergence since the 
numerical scheme being used is of second order, its invariantization was shown to preserve this 
order and the quadratic interpolation and its invariantization is also of second order. 

8.1 Invariant scheme without re-mapping 

In this test run we use the scheme (10) without re- mapping. As a result, the solution is evolving 
along the trajectories of the grid equation (11). Since the spacing between trajectories is not 
constant, i.e. h + and h~ are changing in time, we choose to plot the error versus 1/N. In Fig. 2, 
we observe the second order convergence expected. 




Figure 2. Convergence plot for the no-remapping scheme. 

8.2 Invariant scheme with non-invariant quadratic interpolation 

In this scheme we interpolate the solution of the invariant scheme (10) at every step back onto 
the regular grid using standard Lagrange quadratic interpolation (20). As a result, we have 
the solution on a regular grid with step-size h = 1/N. In Fig. 3, we observe the second order 
convergence expected. 
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Figure 3. Convergence plot for the non-invariant quadratic interpolation scheme. 

8.3 Invariant scheme with invariant quadratic interpolation 

In this scheme, we interpolate the solution of the invariant scheme (10) at every step back 
onto the regular grid using the invariant Lagrange quadratic interpolation (21) described in the 
previous section. As for the case above, we have the solution on a regular grid with step-size 
h = 1/N. In Fig. 4, we observe the second order convergence expected. 




h 



Figure 4. Convergence plot for the invariant quadratic interpolation scheme. 

8.4 Linearity preservation in the invariant numerical scheme 

Linearity is not explicitly preserved by the schemes proposed. For this reason it is instructive 
to check whether or not the linearity is preserved in the fully invariant scheme. 
Consider, 

Ut = U X x 
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u (x, t = 0) = (sin (x - 1) + 2) + (cos (x) + 2) , 

which solution we call u exac t. We then solve numerically the following two equations 

11% = u xx with u a (x, 0) = sin (x — 1) + 2, 
u t = u xx with u fc (x, 0) = cos (x) + 2, 

and define u s = u a + u b . 

Fig. 5 depicts the L oo error between u eX act 

and u s . We observe a convergence rate of second 
order. In other words, despite the fully invariant numerical scheme does not explicitly preserve 
symmetry associated with the linear superposition principle, we observe that linearity property 
is preserved approximately to the order of the method. 




h 



Figure 5. Convergence plot for the linearity test using the invariant scheme with invariant quadratic interpolation. 



9 Conclusions 

In this paper we construct invariant discretization schemes using the method of invariantization 
via equivariant moving frames. The advantage of this technique is that it allows one to start with 
a given non-invariant scheme and convert this initial scheme into a finite difference approximation 
of a system of differential equations C that is invariant under the same maximal Lie invariance 
group G (or a suitably chosen subgroup of G) as admitted by C. The possibility of converting 
non-invariant numerical schemes into invariant discretizations may lead to the overly optimistic 
speculation that the schemes constructed by invariantization could be included easily in the 
existing numerical models written for the original scheme. The hurdle preventing this in practice 
is that preserving symmetry groups of systems of evolution equations more complicated than 
scalings or translations requires the usage of moving grids. Converting numerical models that 
use standard discretization schemes based on fixed lattices to (invariant) discretization schemes 
on moving meshes is not an easy task. At the same time, rewriting numerical models from 
scratch for the simulation of involved physical processes using symmetry-preserving schemes 
might be a time-consuming and costly task too and it not certain that this is feasible at all. 
Moreover, it is as of now unclear whether the preservation of symmetries by numerical schemes 
for multi-dimensional systems of partial differential equations gives that much of an added value 
compared to standard schemes that one might justify such an undertaking in practice. 
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This is why one relies on finding methods allowing one to efficiently include invariant dis- 
cretization schemes into existing numerical models without the need to rewrite new models 
from scratch that incorporate the invariance methodology. The method proposed in this arti- 
cle solves this problem by breaking the integration procedure into two steps, the time-stepping 
using the invariant numerical scheme with an invariant numerical grid equation and the pro- 
jection (interpolation) of the results obtained at intermediate grid points to the regular mesh. 
This interpolation can be done in an invariant way by applying the moving frame map used to 
invariantize the initial discretization scheme also to a particular interpolation method. An al- 
ternative is to assemble the invariant interpolation method using the difference invariants which 
follow from applying the invariantization map to the single points of the lattice. Either way, it 
is worthwhile pointing out that interpolations requiring only data given at a single time level 
are already invariant under most symmetry groups as admitted by physical systems of differ- 
ential equations. Thus, invariantization of interpolation formulas will often only lead to minor 
modifications of the initial interpolation method chosen and the influence on the numerical so- 
lution might be rather small. In the numerical tests carried out above for the heat equation, 
the difference in the convergence properties we found when using invariant or non-invariant 
interpolation methods is indeed small although using the invariant interpolation gave slightly 
better numerical results. This is encouraging and the reason why we plan to further investigate 
invariant numerical schemes using the re-mapping procedure. 

We illustrate the evolution-re-mapping strategy by integrating the one-dimensional linear 
heat equation with an invariant numerical scheme. The heat equation has been studied quite 
extensively in light of its invariance properties and in particular it is a standard model for 
the construction of invariant numerical schemes [1, 13, 34]. At the same time, a comprehen- 
sive numerical analysis of such schemes was not given before and thus seems relevant to be 
reported. This is another aim of the present paper. Again, the analysis of numerical properties 
of discretization schemes is considerably easier if one can use non-evolving meshes. 

Further work we intend to do is to employ the evolution-re-mapping strategy to multi- 
dimensional systems of differential equations using both higher-order discretization and inter- 
polation schemes. 
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